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Abstract 


A time dependence of three-dimensional simulation including water phase change and heat transfer of a PEMFC model has been studied. 
The overshoot behavior has been observed during a change in the electrical load during operation with fixed flow rates of hydrogen and air. The 
simulation of 25-cm° active area with a serpentine flow path shows the interactions of the anode and cathode flow streams, the flow through the 
gas diffusion media, and the movement of water through the MEA by electroosmotic and back diffusion forces. The simulation used a commercial 
computational fluid dynamics (CFD) solver, STAR-CD with and add-on PEMFC module, es-pemfc. The operating conditions corresponded to 
101 kPa, 70°C cell temperature, anode and cathode dew points and stoichiometries of 78 °C and 72°C and 1.2 and 2.0 at an initial operating voltage 


of 0.7 V and current density of 0.57 A cm™°?. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Proton exchange membrane fuel cells (PEMFCs) produce 
power for electric drive motors with less pollution. A full-scale 
three-dimensional solution to the time dependent Navier-Stokes 
equations for the flow channel and diffusion layers has been 
developed by including water phase change model [1,2] to 
improve the investigation of transient behavior of the fuel cell 
and its performance, which can be changed to meet the require- 
ments of the amplitude and frequency of the load changes. 

Many researchers have studied the problem of water manage- 
ment and distributions inside the PEM fuel cell especially under 
steady state operation as discussed in Wang [3] and Baschuk and 
Li [4]. Moreover, several groups [5—17] are improving their com- 
putational models to be more realistic, faster in computing, or to 
be able to be used for design improvement of PEMFC. Meng and 
Wang [8] improved their three-dimensional (3D) computational 
fluid dynamics (CFD) model to be more accurate to investigate 
two-phase behavior under different gas utilizations although the 
energy transport was ignored. Kulikovsky et al. [10] presented 
a simplified analytical model to obtain physical parameters and 
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then input those parameters into a quasi-3D model to speed up 
the PEMFC calculation. Oosthuizen et al. [13] used a 3D model 
to study the gas crossover between side-by-side channels under 
different flow rates, channel path length flow fields, and GDL 
porosities. Their results might be useful for flow field design 
improvements of the PEMFC. The multidimensional model cal- 
culations provide details inside the fuel cell on a local level and 
describe distributions of current, heat, and water. Thus, mod- 
eling will help in gaining an understanding of the mechanisms 
inside the fuel cell, aid in experimental data analysis, and identify 
limiting parameters. 

The above mentioned modeling studies simulated a steady 
state operation of the fuel cell but a dynamic simulation is needed 
for PEMFC while operating conditions are changing, e.g. when 
the power required is changing under automotive application. 
Kim et al. [18—20] presented transient experimental results of 
25-cm? active area PEMFC while the cell potential was varied 
with different rate. They also studied the effect of flow field 
design and fuel dilution on dynamic response of PEMFC. Their 
results indicate current overshoot behavior while the system 
requires more power rapidly and undershoot behavior while it 
requires less power. The peak of overshoot and undershoot cur- 
rent density depend on flow filed design, range of power change, 
gas component, and manifold configuration. Many experimen- 
talists continue working on PEMFC research under dynamic 
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Nomenclature 


br 


activity of water in stream E 

specific surface area of the control volume (c.v.) 
(m!) 

surface area of control volume (m7) 
condensation rate (s7!) 

concentration of water vapor at kth interface of 
the membrane (mol m~?) 

concentration of water liquid at kth interface of 
the membrane (mol m7?) 

hydraulic diameter of flow channel (m) 

diffusion coefficient of Hz in liquid water film 
(6.3E—9 m? s7!) 

binary diffusion coefficient of species n in mixture 
jm? s~!) 

diffusion coefficient of O2 in liquid water film 
(2.4E—9 m? s7!) 

diffusion coefficient of water (m? s7!) 

Faraday constant (96,487 C mole-of-electrons™ !) 
enthalpy of vaporization for water (kJ kg!) 


hformation enthalpy of water formation (kJ kmol~!) 


Ay, 1 


Ho), 


Henry’s law constant for hydrogen in liquid water 
film (8.9E—9 Pa) 

Henry’s law constant for oxygen in liquid water 
film (2.12E—10 Pa) 

local current density (A m~?) 

exchange current density (A m~?) 

mass of species n (kg) 

mass fraction of the species n in stream k 

mass fraction of liquid water 

equivalent weight of a dry membrane (kg mol!) 
molecular weight of species n (kg mol!) 
electro-osmotic drag coefficient (number of water 
molecules carried per proton) 

pressure (Pa) 

partial pressure of species (Pa) 

vapor pressure of water in stream k (Pa) 

volume flow rate (m? s7!) 

universal gas constant (8.314 J mol`! K7!) 
source term 

time (s) 

liquid water film on the anode (m) 

liquid water film on the cathode (m) 

membrane thickness (m) 

temperature (K) 

velocities in x, y, and z directions, respectively 
(ms~') 

cell open-circuit voltage (V) 

cell voltage (V) 

channel length measured from anode inlet (m) 
mole fraction of species n 

mole fraction of water in stream k 


Greek symbols 
a net water flux per proton flux 
Oy anode transfer coefficient 


Oe cathode transfer coefficient 

Be permeability in the E direction 

E porosity of gas diffusion layer 

n over potential for oxygen reaction (V) 
A water content in the membrane 

u dynamic viscosity (kg s~! m7?) 

p density of the mixture (kg m7?) 
density of a dry membrane (kg m~?) 
Om membrane conductivity (Q-! m7!) 


Subscripts and superscripts 


a anode 

CG cathode 

e electrochemical reaction 
H: hydrogen 

K anode or cathode 

l liquid 

N2 nitrogen 

O2 oxygen 

p phase change 

sat saturated 

v vapor 

W water 

E dummy variable for direction x, y, or z 


response including stack design, contaminants, flow-field, etc. 
For example, Hamelin et al. [21] revealed three different tran- 
sient responses on Ballard fuel cell stack model MK5-E, which 
has maximum output of 10kW. Amphlett et al. [22] gave the 
experimental results of a hybrid PEMFC/battery system. They 
showed how their hybrid system and components interacted dur- 
ing charging or discharging. Mohtadi et al. [23] showed dynamic 
response of current density when PEMFC was exposed by H2S 
while analyzing the effect of temperature on the adsorption rate 
of H2S on Pt anodes. Emonts et al. [24] studied the dynam- 
ics of a PEMFC, a compact methanol reformer, and the choice 
of a short-term storage system for an automotive application. 
Only few research groups have included transient simulation of 
PEMECs in their studies. For instance, Wang and Wang [25] pre- 
sented dynamics study on PEMFC using commercial CFD code. 
Their geometry was single straight channel flow-field that is not 
normally used in fuel cell application. Therefore, the effect of 
pressure driven flow underneath the gas diffusion layers has been 
ignored. Shimpalee et al. [26,27] showed numerical solutions 
of transient response of 10-cm? serpentine flow-field PEMFC 
under both excess and starve fuel/air. They concluded that the 
peak of overshoot depends on the rate of voltage change and 
utilization of reactants. The overshoot behavior is a result of 
local non-uniformity in distributions. However, their model did 
not take into account of temperature effect and phase change of 
water. Therefore, their results might be over predicted. 
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Current collector 
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Fig. 1. The picture shows schematic of fuel cell assemble and actual flow-field plate with gas channel of 25-cm? PEMFC. There are 30 straight channels connected 
in a triple pass serpentine fashion. Anode side and cathode side flow channels are symmetric and placed property aligned (non-staggered) on top of each other. 


In this study, the three-dimensional model of PEMFC with 
25-cm? reacting area as shown in Fig. 1 was developed with 
the inclusion of time dependent analysis. The transient model of 
Shimpalee et al. [26,27] was used with the extension of energy 
and water phase change equations. Different rates of voltage 
change were chosen to study their effects on the PEM fuel cell 
performance at each time interval. Moreover, the detail of local 
current density, water distribution, and reactant concentration 
are investigated and compared with previous works [26,27]. 


2. Numerical procedure 


Fig. 2 shows model geometry that created from Fig. 1. It 
has triple pass channel with 11 serpentine fashions. The geom- 
etry of fuel cell arrangement modeled in this work consists of 
two flow channels separated by diffusion layers and the mem- 
brane electrode assembly (MEA). Between diffusion layers and 
MEA, there are sub-gaskets. Therefore the active area of MEA 
is reduced from 25-cm? to 20-cm? as shown in Fig. 2b. There are 
triple pass channel with eleven serpentine passes in the flow path, 
so that the flow is approximately 55 cm long in the axial direc- 
tion with 0.1 cm (height) x 0.08 cm (width) cross-section flow 
area in each channel. Each diffusion layer has dimensions of 
0.025 cm (height) x 5.00 cm (width) x 5.00 cm (length). A total 
of 734,940 cells (elements) were used to model the fuel cell. 

The equations used in this work are shown in Tables 1 and 2. 
These time dependent equations have been extended from Shim- 
palee et al. [26,27] to include the energy equation, and a water 
phase change model, where homogeneous two phase flow was 
assumed. In the water phase change model, when the local activ- 
ity of water exceeded 1.0, water vapor was condensed to form 
liquid water until the local activity equaled 1.0. Conversely, if 
liquid water was present and the local activity of water dropped 
below 1.0, then liquid water was evaporated until the local activ- 
ity equaled 1.0. When liquid water condensed in the region 
adjacent to the electrode surface, it is assumed to form a liq- 
uid film on the electrode surface. In the regions where there 
was a liquid film on the electrode, hydrogen and oxygen were 
required to dissolve in the liquid film and diffuse through the 
film to the electrode surface in order to react. Henry’s law was 
used to calculate the solubility of the gases in the liquid. The 
thickness of this liquid film depends on the rate of condensa- 


tion/evaporation and the production of water by electrochemical 
reaction as discussed in the appendix of Lee et al. [5] and Shim- 
palee and Van Zee [2]. Note that since the solubility of O2 in 
the liquid film is low, the solution procedure includes transient 
calculation (see Xo, (x, y, z, t) in Eq. (30)) in the water film on 
the catalyst surface. This transient concentration and water film 
thickness corresponds to a lost in available surface reaction area 
when the film is thick or three phase boundary. 

A control volume technique based on a commercial flow 
solver, STAR-CD 3.26, was used to solve the coupled governing 
equations [28] with transient analysis. This software was used 
with an add-on tool called Expert System for PEMFCs (es-pemfc 
version 2.20) that provided the source terms for the species trans- 
port equations, the phase change equations for water, and the heat 


(b) Graphite Block 


Sub-gasket 


Flow channel covered 
by Sub-gasket (Dotted line) 


Flow channel uncovered 
(Solid line) 


Fig. 2. (a) The geometrical model of the complete 25-cm? fuel cell shown with- 
out graphite current collector. There are 33 straight channels connected in a 
triple pass serpentine fashion shown in Fig. 1. (b) Schematics of sub-gasket and 
flow channel (20-cm? active area due to sub-gasket). 


S. Shimpalee et al. / Journal of Power Sources 167 (2007) 130-138 


Table 1 
(a) Governing equations and (b) source terms for governing equation in (a) 
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Governing equations Mathematical expressions 


(a) 
d a a a 
Conservation of mass a t (ou) t CH H (ow) = Sm 2) 
dr ox dy az 
dou a( pu) d(pu) (pu) oP d ou d ou 
u v Lu = d H: H F Spx 
o 3 ax e dy à az Pi oe a ig ae 
Momentum transport Pr u (ov) v (ov) + Ww (ov) = H- Ga Gë d dÉ z) + Spy (3) 
ðt ox dy d du ox ax H May az dz ` 
dou d(pw) dou) dou) ðP z(a? ow 7 E (+; 2) d ðw ) 
Lu Lu Hw = H H- + Spz 
at ax y z Z a dy az d 
d d d d 7 Jx m A d Jz 
Hydrogen transport (anode side) A u "ri v a w "ni L Lan i (ym) ( et) SH (4) 
D y z "` 
d ; d d a d Z d b j d b 
Water vapor transport PM wy i (pmwy) o (emwy) ‘ai (omwy) 8 ( cal ( oO (Je, a EECH (5) 
at ax dy dz ax dy R 
d a d d a(J, a(Jy a(Jz 
Water liquid transport PM w] e (omw) Ly (omw) L (omw) _ ( x, wi) ( yw) ( 5 Sn (6) 
a S „2 n A a a i ad ar d d 
TE ats) pmo, (emos) Loo: M (pmo) _ ( v.02) e Lo) So, D) 
ot dy dz ox dy d S 
È 8 dph , oul) | wer, d(pwh) d (i =) , 9 (1 =) E (1 =) LK A8 (8) 
nergy equation ap ae By z Ee ax) T By y) z za hp + She 
Non-zero volumetric source terms 
(b) 
Conservation oi mass Sm = SH, + Swvp + Swip + Swve at anode GDL/MEA interface 9) 
Sm = So, + Swvp + Swip t Swe at cathode GDL/MEA interface 
Momentum transport Spx = E Spy = , Sp: = = at GDLs (10) 
Bs i z 
yt 3 
Hydrogen transport (anode side) Sm, = — “ > M H Acy at anode GDL/MEA interface a1) 
M) mass, vz, TI Muotn psat T.t)— t 
Water vapor transport Be, = ER otu paral eae vB 1) — Pov y zt) * rat entire fluid domain (12) 
E 1 — (PACT, SSC ,Z,t)) P(x, y, z, t) 
a(x, y, t) 
Sawve = — ——— Pa , t)Mp,0 Acv at anode GDL/MEA interface (13) 
1+ 2a(x, y. 
Scwve = FOO I, y, 1)My,0 Acy at cathode GDL/MEA interface (14) 
Water liquid transport Swip = —Swvp at entire fluid domain (15) 
I t 
Oxygen transport (cathode side) So. = — “ = ) Mo, Acy at cathode GDL/MEA interface (16) 
Energy equation Heat source by energy losses: 
1 
She = Mformation * l= I(x, y, t) Ae) — U(x, y, t)Vcett Acy) at cathode GDL/MEA interface (17) 
Heat source by phase change: 
Shp = Swp * hfg at entire fluid domain (18) 


generation equations [2,29]. Also, es-pemfc accounted for the 
flux of protons and water across the membrane [2,29]. The mate- 
rial properties and operating conditions of this work are shown 
in Table 3. The transient load profile consists of changing the 
cell voltage from 0.7 V at steady state to 0.5 V with the rate of 
1.0 V s7! as shown in Table 4. 

The transient solution procedure used in the STAR-CD flow 
solver is based on a PISO algorithm [28] where a predictor- 
corrector strategy is used. For each time step, a predictor, 
followed by a number of correctors, iteratively solves the lin- 
ear equation sets for each main dependent variable. The mixture 
properties at each control volume are calculated based on the 
local species concentrations. The anode side gas mixture con- 
tains hydrogen and water. And, the cathode side gas mixture 
contains oxygen, water, and nitrogen. Therefore, the density 
and viscosity of the two flow channels are different and they 
vary from at each location and the newly calculated values of 
the variables prevail throughout the time step At. 


3. Results and discussion 


Fig. 3 shows how the averaged current density changes 
when the cell voltages are changed from 0.7 V to 0.5 V using 
transient load profile from Table 2. This figure indicates over- 
shoot behavior. The averaged current density increases from 
0.57 A cm7? to 1.23 A cm~? when the cell voltage is 0.52 V at 
time = 0.23 s and then the current density slowly decreases and 
reaches 1.15 A cm7? at time=0.35s. The height of the over- 
shoot can be changed by the rate of cell voltage change and by 
the flow rates as already discussed by our previous work [26] for 
the case of excess stoichiometry but here we fix the rate of volt- 
age change and flow rate and study the dependent variables in 
order to understand their interactions during current overshoot 
behavior. 

Fig. 4 shows the local current density distributions on the 
membrane surface at four different time steps and cell voltages. 
At the steady state current density of 0.57 A cm~? correspond- 
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Table 2 
Constitutive equation for modeling electrochemical effects 
d 
Diffusion mass flux of species n Jen = —pDeEn SE (19) 
i KE 0g 
in é direction 
F CwelX, y, t) — Cwalx, y, t 
Net water transfer coefficient per a(x, y, t) = nat, y, t) — ———— Dy (x, y, D EE EE (20) 
I(x, y, t) tm 
proton 
Water content in the membrane à = 0.043 + 17.8la, 39.852 d 36.043; 0 < da < 1 = 14 + 1.4(aa — 1); 1 < aa < 3 (21) 
Electro-osmotic drag coefficient na(x, y, t) = 0.002914? + 0.05 — 3.4 x 107! (22) 
1 
Water diffusion coefficient Datz, y, t) = Dy exp (2416 E m =)) D = 107!?, à < 2; D} = 107191 + 204 — 2), 2 < A < 3; D} = 
x,y, 
107103 — 1.67(A — 3), 3 < à < 4.5; Dy = 1.25 x 10719, à > 4.5 (23) 
Water vapor concentration for Cwk(x, y, t) = ee (0.043 + 17.8aK 39.84% d 36.0a%); aK aS Pindiy (14 + 1.4(ax — 1)); for ax > 
anode and cathode surfaces of m,dry "dn 
the MEA 1, where K = aorc (24) 
X 2, D) P(x, y, z, t 
Water activity aK = SE? x A (25) 
PakO D 
a E 
Local current density I(x, y,t) = m y Im. Maa — 9.1, (26) 
m 
GES Ma dr 1 1 2 
Local membrane conductivity Om(x, y, t) = | 0.00514—— Cwalx, y, t) — 0.00326 | exp | 1268 | — — ———— x 10 (27) 
Pm, dry 303 T(x, Y; t) 
i RT(x, y, cathode -GDL/MEA, t) I(x, y, t)P(x, y, cathode -GDL/MEA, t) 
Local over potential n(x, y, t) = In 
oP Ins, Po, (x, y, catbode GDL /MEA. t) 
RT(x, y, cathode.GDL/MEA, t) i I(x, y, t)P(x, y, anode- GDL/MEA, t) 28) 
Ol, F "| Ton, Pigs y, anode_GDL/MEA, 1) 
m mass, 
Water film thickness tk = at, massa) (29) 
EPwiATeacy 
Gas solubility Po, (x, y, z, t) = Xo, x, y, z, OPC, y, z, t) at cathode GDL/MEA interface (30) 
I(x, y, t) Xo, (x, y, z — te) PC%, y, z — thes Don — Xor, y, z, f) 
———— Mo, = pmo, Do, - where z = 
4F thc 
cathode GDL/MEA interface (31) 
Puy (x, y, z, t) = Xp, (x, y, z, PC, y, z, t) at anode GDL/MEA interface (32) 
EE Xm (xX, Y, Z + tia, DPG, y, z + tha, Nol — Ent, y, z, £) 
e Mn, = pmp, Dy, 1 ` where z = 
2F 5 tha 
anode GDL/MEA interface (33) 


Table 3 
Material properties and operating condition at 0.7 V corresponding to stoichiom- 
etry of 1.2/2.0 at 0.57 A cm~? 


Anode channel inlet Velocity (ms~!) 2.23 
conditions Mole fraction of H? 0.53 
Mole fraction of H2O 0.47 
Dew-point temperature (°C) 78 
Cathode Velocity (ms~!) 6.8 
chan- Mole fraction of O2 0.14 
nel Mole fraction of N2 0.50 
inlet Mole fraction of H2O 0.36 
con- Dew-point temperature (°C) 72 
Operating Cell temperature (°C) 70 
con- Operating pressure (kPa) 101 


Den Collector Thermal conductivity (W m7! K7!) 5.7 


GDLs Thickness after compressed (um) 250 
Permeability (m?) 3.3E—15 
Porosity after compressed (%) 70 
Diffusion adjustment (%) 50 
Thermal conductivity (W m7! K7!) 0.21 
Membrane Gore PRIMEA®5561 


ing to the cell voltage of 0.7 V and t=0.10s, the local current 
density is higher at the inlet region and decreasing along the flow 
channel towards the outlet due to the consumption of reactants 
and the reduction of anode water activity. The maximum and 
minimum local densities are ~0.60 A cm7? and ~0.50 Acm72, 


1.3 0.75 
Overshoot —>, 
1.24 
+0.70 
E 1.14 
+0.65 
z 1.04 S 
D 
SEH 0.60 & 
S & 
Cosi S 
c -0.55 = 
D E 
3 0.74 O 
+ 0.50 
0.64 
0.54 : : voltage decreasing rate TOV : f Lo.45 
0.00 0.0 0.10 0.145 0.20 025 0.30 0.35 0.40 
Time (s) 


Fig. 3. Average current density for change from 0.7 V to 0.5 V. 
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t=0.30s, Veen = 0.5V, lava = 1.19 A/cm? 


t=0.35s, Veen = 0.5V, lave 


2.180 
2.057 
1.934 
1.811 
1.689 
1.566 
1.443 
1.320 
1.197 
1.074 
0.9514 
0.8286 
0.7057 
0.5829 
0.4600 


=1.15 A/cm? 0,405, Veet = 0.5V, Iavg = 1.13 Alem? 


Fig. 4. Local transient current density (A cm~*) contours on MEA surface at different times and cell voltages. 


respectively, which represents small variation over the entire 
active cell area. When the cell voltage is dropped to 0.52 V 
at time = 0.23 s, the local current density significantly increases 
over the membrane surface however, the contour pattern is simi- 
lar to when cell voltage is 0.7 V. The highest local current density 
is ~2.18 Acm~? and the lowest value is ~0.95 A cm". When 
the time increases to 0.30 and the cell voltage is 0.5 V, there is 
more non-uniformity of local current density distribution. The 
local current reduces considerably on the MEA region from 
the center of membrane surface toward the outlet due to strong 
decreasing of oxygen concentration. The highest local current 
density is located around the inlet with a value of 2.0 A cm~? 
and the lowest local current density is at the outlet with a value of 
0.58 Acm~*. When PEMFC performance reaches steady state 
at cell voltage of 0.5 V and time is greater than 0.35 V, the 
decreasing of local current density from the center region of 
the membrane surface towards the outlet region is more signifi- 
cant. The current density in this case varies from 2.0 A cm7? to 
0.46 A cm". which represents considerable variations along the 
membrane surface. For the case under consideration, at steady 
state cell voltage of 0.5 V, the local current density does not sig- 
nificantly vary over the first half of the membrane surface, but 
it sharply decreases over the second half of membrane surface. 
This is because the high reaction rate regions (i.e., first half of the 
membrane surface) accompany excessive hydrogen and oxygen 
consumptions and therefore they depleted in these regions (i.e., 
second half of the membrane surface). 

To further describe the overshoot phenomena, the profile plots 
across MEA width and serpentine flow path at x=2.5cm as 
shown in Fig. 5 are examined. Fig. 6 shows the transient profiles 
of local oxygen mole fraction across the MEA width and channel 


as shown in Fig. 5. At cell voltage of 0.7 V and steady state, the 
oxygen mole fraction decreases along the electrode width from 
inlet toward the outlet due to its consumption from electrochem- 
ical reaction, which is consistent with the local current density 
distribution from Fig. 4. When cell voltage is reducing to 0.58 V 
at time =0.2s, its distribution shows more non-uniform due to 
the higher reduction of oxygen concentration by higher reaction 
rate. The much more non-uniform distribution of oxygen mole 
fraction is revealed at cell voltage of 0.52 V at time =0.23 s, 
especially in the first half of the cell where the reaction rate 


5.0 cm > 


i Outlet 


Section Cut 


j Inlet 


Fig. 5. Top view projection of serpentine flow-field showing location of analysis. 
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Table 4 

Cell voltage changes from 0.7 V to 0.5 V in 0.2s 

Time (s) 0.10 0.125 0.15 0.17 
Cell voltage (V) 0.70 0.695 0.685 0.67 


Predicted current density (A cm7?) 0.56 0.58 0.62 0.69 


0.18 


0.19 0.20 0.21 0.23 0.24 0.26 0.28 0.30 
0.63 0.60 0.57 0.52 0.51 0.505 0.5025 0.50 
0.76 0.85 0.98 1.11 1.23 1.22 1.21 1.18 1.15 
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E 
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D 
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l hap tere 
0.024 Z MIN = = 
t=0.36s, veel um" 1 Une d tege Keng 
VV 
0.00 + T T T: T p mf T T i 
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Y-axis (cm) 


Fig. 6. Variation of oxygen mole fraction along channel width at different times 
and cell voltages. 


is higher and it has much lower concentration in the second 
half of the cell. Again at this time, the PEMFC gives the maxi- 
mum average current density. When the cell voltage is reaching 
0.5 V after time greater than 0.23 s, the oxygen concentration 
is depleting in the second half of the cell Oe, a mole fraction 
between 0.02 and 0.005) and this limits the current density in 
these locations. Therefore, the current density becomes decreas- 
ing thus indicating overshoot behavior. When cell voltage of 
0.5 V reaches steady state performance, the slope of oxygen 
mole fraction profile becomes negatively steeper in the first half 
of the cell and remains minimal value in the second half of the 
cell. In summary, when the cell voltage is decreasing, the rate 
of gas consumption is increasing and the oxygen concentration 
decreases. 
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Fig. 7. Variation of cathode over potential (V) along channel width at different 
times and cell voltages. 


Fig. 7 presents transient response of cathode over potential 
or cathode polarization that is a function of partial pressure of 
oxygen. As already indicated in Fig. 6 as oxygen mole fraction 
response with time and cell voltage change, those profiles affect 
cathode kinetic significantly. At time =0.1 s and cell voltage is 
0.7 V, there is not significant polarization at cathode surface 
because there is plenty of oxygen available in the PEMFC as 
shown in Fig. 6. When current is increasing and reaches the 
peak overshoot (cell voltage = 0.52 V and time is 0.23 s), cath- 
ode over potential or polarization increases with about the same 
slope or rate along electrode width (Y-axis) as the initial point of 
0.7 V. This is because the consumption of oxygen is increasing 
but its concentration has not depleted. After t= 0.23 s and the cell 
voltage is remaining constant at 0.5 V with constant flow rate of 
reactants, the cathode polarization is increasing especially at Y- 
axis is greater than 2.5 cm (second half of MEA surface), which 
is consistent with the depletion of oxygen concentration as dis- 
cusses in Fig. 6. This profile reaches steady state of Veeg =0.5 V 
at t=0.35 s. This will drop performance or current density thus 
indicating overshoot phenomenon. 

Fig. 8 shows transient response of liquid water film thick- 
ness profiles along the MEA width at the cathode side. This 
figure shows that the liquid water film is increasing from 
inlet toward outlet due to the production of water vapor by 
electrochemical reaction and then water vapor can be con- 
densed to liquid water if water activity is greater than 1.0. 
This figure also shows the accumulation of liquid water from 
time =0.10s at 0.7 V to time=0.35 at 0.5 V. This accumula- 
tion of liquid water could decrease the PEMFC performance 
by creating the high gases resistance known as water flood- 
ing. However in this condition, the liquid water presented in 
this system seems to be insignificantly to decrease the PEMFC 
performance. 
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Fig. 8. Variation of liquid water film thickness (micron) along the centerline on 
cathode electrode at different times and cell voltages. 
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Fig. 9. Average cell current, cathode over potential, ohmic loss, and temperature 
response during a transient load change. 


Fig. 9 presents the transient average temperature on the anode 
MEA surface, cathode over potential, and ohmic loss correspond 
to particular average current density and the cell voltage at each 
time step. This results shows that temperature increases when 
the cell voltage is decreased. This is because the reduction of cell 
voltage escalates the reaction rate thus increasing heat genera- 
tion caused by electrochemical reaction. Note that the increase in 
temperature results in lower membrane water content, thus drop- 
ping in ionic conductivity. Further, the rate of temperature rise 
is much slower than the rate of current density change. There- 
fore, the overshoot in temperature is not as significant as the 
current density overshoot in this study. For the cathode over 
potential profile, it increases when the current density is increas- 
ing, which is similar to the ohmic potential profile. However, 
after the current density reaches the peak of overshoot, cathode 
over potential continues increasing due to the loss of oxygen 
concentration near the exit region as explained earlier. Conse- 
quently, current density starts to decrease, thus causing the lower 
temperature. Meanwhile, the ohmic loss is slightly decreasing 
due to the increasing in membrane water content because of the 
lower temperature. Therefore, the contribution of ohmic loss or 
membrane conductivity and temperature in this particular oper- 
ating condition could be the controlling factor in the peak and 
magnitude of overshoot. 

Fig. 10 presents the comparison of transient behavior between 
the transient current results of single phase isothermal model 
[27] on 10-cm? active area and phase change model on 25-cm? 
active area. For better comparison, the current densities shown in 
this figure need to be dimensionless with average current density 
of each time step (/(f)) is divided by average steady state current 
density at cell voltage of 0.5 V Usteady state at 0.5V) as equation 
presented below. 

W= — a) 
I steady state at 0.5 V 

The result reveals that the current density of single phase, isother- 

mal model has much higher overshoot than the result of water 

phase change, non-isothermal model. Moreover the rate of cur- 

rent density reduction after overshoot point of single phase, 

isothermal model is faster than phase change model. This could 
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Fig. 10. Comparison of transient behaviors between single phase, isothermal 
model and water phase change, non-isothermal model. 


be the effect of temperature and water phase change that included 
into this recent model. 


4. Conclusions 


A three-dimensional including water phase change and 
energy balance with transient simulation of PEM fuel cell was 
studied at fixed flow rates of feed gases. The flow rates cor- 
respond to normal fuel and oxygen utilization (1.2/2.0 stoich) 
at the high cell voltage and a minimal operating condition at 
the low cell voltage. The cell voltage was changed from 0.7 V 
to 0.5 V at a somewhat rapid rate of 1.0 V s7! and overshoot 
behavior in the current was observed and explained in terms of 
the distributions of current density, oxygen concentration, liquid 
water, and temperature. 

When the cell voltage is decreased, the current density 
increases and reaches the peak of overshoot while the capacity of 
the reactants in the flow channels is used. This steady state values 
correspond to a highly non-uniform local current density distri- 
bution, dependent on local concentrations of the reacting gases. 
For this condition, oxygen appears to control the distribution. 
Single phase and isothermal model shows higher overshoot than 
water phase change and non-isothermal model because temper- 
ature rise affects the peak of overshoot. Therefore, water phase 
change and non-isothermal effects are important in a transient 
PEMFC model. Moreover, the controlling mechanism of over- 
shoot is dependent on the gases utilizations and the rate of load 
change. These results may help effort in the design of the flow- 
field and manifold configuration. These design changes may 
improve stack performance. 
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